{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d405c15d-533e-4a3c-a9a9-25792f41a633",
   "metadata": {},
   "source": [
    "### Enrichment of misexpression-associated SVs in regulatory annotations \n",
    "\n",
    "* 3D genome architecture \n",
    "    * TAD boundaries from: \n",
    "        * GM12878 (shared across multiple cell types) \n",
    "    * A/B compartments \n",
    "        * GM12878 cell line \n",
    "* CTCF-binding sites \n",
    "    * All CTCF-only \n",
    "    * B-cells \n",
    "    * Neutrophils \n",
    "    * CD14 monocytes \n",
    "* Chromatin features \n",
    "    * Chrom HMM \n",
    "* CpG islands "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b70790a6-1632-48e0-a9db-2207a620ddec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np\n",
    "from pybedtools import BedTool\n",
    "from io import StringIO\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.discrete import discrete_model\n",
    "import sys\n",
    "from collections import Counter\n",
    "from scipy.stats import fisher_exact\n",
    "from patsy import dmatrices\n",
    "from pathlib import Path\n",
    "import math\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.style as style\n",
    "from matplotlib import pyplot\n",
    "from functools import reduce\n",
    "from statsmodels.stats import multitest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "44ee9daa",
   "metadata": {},
   "outputs": [],
   "source": [
    "wkdir = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/\"\n",
    "wkdir_path = Path(wkdir)\n",
    "# inputs \n",
    "sv_info_path = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/info_table/final_sites_critical_info_allele.txt\"\n",
    "all_vrnts_path = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/vrnt_id_in_window_cntrl_misexp_genes.txt\")\n",
    "all_vrnts_bed_path = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/vrnt_id_in_windows_misexp_genes.bed\")\n",
    "misexp_vrnts_path = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/vrnt_id_misexp_tpm_zscore_median.txt\")\n",
    "# shared TAD boundaries across cell lines \n",
    "shared_tad_boundaries_path = wkdir_path.joinpath(\"reference/4d_nucleome/shared_boundaries/4DNFIVK5JOFU_imr90_huvec_hnek_hmec.bed\")\n",
    "# ENCODE cCREs \n",
    "encode_c_cres_dir = wkdir_path.joinpath(\"reference/encode/encode_c_cres\")\n",
    "# CpG islands \n",
    "cpg_isl_bed_path = wkdir_path.joinpath(\"reference/cpg_islands/cpgIslandExt.bed\")\n",
    "# PBMCs chromHMM \n",
    "pbmcs_chromhmm_15marks_path = wkdir_path.joinpath(\"reference/chromhmm/E062_15_coreMarks_hg38lift_mnemonics.bed\")\n",
    "# A/B compartments in GM12878\n",
    "gm12878_ab_path = wkdir_path.joinpath(\"reference/4d_nucleome/gm12878_hi_c/compartments/4DNFILYQ1PAY.bg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9d195f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# output\n",
    "out_dir = wkdir_path.joinpath(\"5_misexp_vrnts/functional\")\n",
    "out_dir_path = Path(out_dir)\n",
    "out_dir_path.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "75eeff15",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of variant IDs in windows: 20255\n",
      "Number of unique misexpression-associated SVs: 105\n"
     ]
    }
   ],
   "source": [
    "# load bed file with variants in windows \n",
    "vrnts_in_windows_bed = BedTool(all_vrnts_bed_path)\n",
    "vrnts_in_windows_bed_sorted = vrnts_in_windows_bed.sort()\n",
    "# all variants in windows df\n",
    "all_vrnts_in_windows_df = pd.read_csv(all_vrnts_path, sep=\"\\t\", header=None).rename(columns={0:\"vrnt_id\"})\n",
    "all_vrnts_in_windows = all_vrnts_in_windows_df.vrnt_id.unique()\n",
    "print(f\"Number of variant IDs in windows: {len(all_vrnts_in_windows)}\")\n",
    "\n",
    "# load misexpression-associated variants \n",
    "misexp_vrnts_ids = pd.read_csv(misexp_vrnts_path, sep=\"\\t\", header=None)[0].astype(str).unique()\n",
    "print(f\"Number of unique misexpression-associated SVs: {len(misexp_vrnts_ids)}\")\n",
    "all_vrnts_in_windows_df[\"misexp_uniq\"] = np.where(all_vrnts_in_windows_df.vrnt_id.isin(misexp_vrnts_ids), 1, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1b4fad00-ad84-4098-a18d-d730d438c58b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gm12878_shared\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "### shared TAD boundaries\n",
    "\n",
    "all_chrom_vrnts_tad_count_dist = []\n",
    "\n",
    "bed_intersect_cols = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\", 4:\"boundary_chrom\", \n",
    "                      5:\"boundary_start\", 6:\"boundary_end\", 7:\"boundary_no\", 8:\"boundary_score\", 9:\"overlap\"}\n",
    "\n",
    "window_intersect_cols = {0:\"boundary_chrom\",1:\"boundary_start\", 2:\"boundary_end\", 3:\"boundary_strength\", 4:\"boundary_score\", \n",
    "                         5:\"vrnt_chrom\", 6:\"vrnt_start\", 7:\"vrnt_end\", 8:\"vrnt_id\"}\n",
    "\n",
    "# get overlap for K562, GM12878 and shared GM12878 boundaries \n",
    "cell_type = \"gm12878_shared\"\n",
    "window = 5000\n",
    "vrnt_tad_features_df = all_vrnts_in_windows_df.copy()\n",
    "print(cell_type)\n",
    "# TAD boundaries bed file \n",
    "tad_boundaries_path = shared_tad_boundaries_path\n",
    "tad_boundaries_bed = BedTool(tad_boundaries_path)\n",
    "input_tad_bed = tad_boundaries_bed.sort()\n",
    "\n",
    "# identify SVs that overlap TAD boundaries exactly\n",
    "sv_intersect_tad_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(input_tad_bed, wo=True)))\n",
    "sv_intersect_tad_df = pd.read_csv(sv_intersect_tad_str, sep=\"\\t\", header=None).rename(columns=bed_intersect_cols).astype({\"vrnt_id\":str})\n",
    "sv_intersect_tad = sv_intersect_tad_df.vrnt_id.unique()\n",
    "\n",
    "# annotate variants intersecting TAD boundaries and window around TAD\n",
    "vrnt_tad_features_df[f\"{cell_type}_intersect_tad_boundary\"] = np.where(vrnt_tad_features_df.vrnt_id.isin(sv_intersect_tad), 1, 0)\n",
    "# check for NaNs \n",
    "print(vrnt_tad_features_df.isnull().values.any())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "017810fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load A/B compartments in GM12878\n",
    "gm12878_ab_bed = BedTool(gm12878_ab_path)\n",
    "gm12878_ab_input = gm12878_ab_bed.sort()\n",
    "\n",
    "# intersect with A/B compartment scores \n",
    "ab_bed_intersect_cols = {0: 'vrnt_chrom', 1: 'vrnt_start', 2: 'vrnt_end', 3: 'vrnt_id', \n",
    "                           4: 'compartment_chrom', 5: 'compartment_start',6: 'compartment_end',\n",
    "                           7: 'score', 8: 'overlap'}\n",
    "sv_intersect_ab_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(gm12878_ab_input, wo=True)))\n",
    "sv_intersect_ab_df = pd.read_csv(sv_intersect_ab_str, sep=\"\\t\", header=None).rename(columns=ab_bed_intersect_cols).astype({\"vrnt_id\":str})\n",
    "# label compartment types \n",
    "conditions = [(sv_intersect_ab_df.score >= 0) & (~sv_intersect_ab_df.score.isnull()), \n",
    "              (sv_intersect_ab_df.score < 0) & (~sv_intersect_ab_df.score.isnull()), \n",
    "              (sv_intersect_ab_df.score.isnull())]\n",
    "values = [\"A\", \"B\", \"Unassigned\"]\n",
    "\n",
    "sv_intersect_ab_df[\"compartment_type\"] = np.select(conditions, values)\n",
    "vrnt_compartment_features_df = all_vrnts_in_windows_df.copy()\n",
    "for overlap in [\"A\", \"B\"]:\n",
    "    vrnt_overlap_compartment = sv_intersect_ab_df[sv_intersect_ab_df.compartment_type == overlap].vrnt_id.unique()\n",
    "    vrnt_compartment_features_df[f\"{overlap}_overlap\"] = np.where(vrnt_compartment_features_df.vrnt_id.isin(vrnt_overlap_compartment), 1, 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b790e13-bf60-4597-a21c-a4a637b11e70",
   "metadata": {},
   "source": [
    "### CTCF elements \n",
    "\n",
    "* cCREs CTCF: all cell-types \n",
    "    - CTCF-only elements from all cell-types \n",
    "* Specific primary cells: \n",
    "    * CD14 monocytes \n",
    "    * B-cells \n",
    "    * Neutrophils "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0dd138ed-2461-4efb-862b-2d7b60339d28",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "all: number of unique variants in dataframe: 20255\n",
      "CD14_monocyte: number of unique variants in dataframe: 20255\n",
      "B_cell: number of unique variants in dataframe: 20255\n",
      "Neutrophil: number of unique variants in dataframe: 20255\n"
     ]
    }
   ],
   "source": [
    "bed_intersect_cols = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\"}\n",
    "encode_c_cres_dir_path = Path(encode_c_cres_dir)\n",
    "\n",
    "bed_path_list = [\"GRCh38-cCREs.CTCF-only.bed\", \n",
    "                 \"ENCFF389PZY_ENCFF587XGD_ENCFF184NWF_ENCFF496PSJ.7group.bed\", \n",
    "                 \"ENCFF035DJL.7group.bed\", \n",
    "                 \"ENCFF685DZI_ENCFF311TAY_ENCFF300LXQ.7group.bed\"\n",
    "                ]\n",
    "\n",
    "cell_types_list = [\"all\", \n",
    "                   \"CD14_monocyte\", \n",
    "                   \"B_cell\", \n",
    "                   \"Neutrophil\"\n",
    "                  ]\n",
    "\n",
    "vrnt_ccre_features_df = all_vrnts_in_windows_df.copy() \n",
    "for cell_type, bed_file in zip(cell_types_list, bed_path_list): \n",
    "    ctcf_peaks_path = encode_c_cres_dir_path.joinpath(f\"{bed_file}\")\n",
    "    peak_bed = BedTool(ctcf_peaks_path).sort()\n",
    "\n",
    "    all_chrom_peak_info = []   \n",
    "    # identify SVs that overlap peaks\n",
    "    sv_intersect_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(peak_bed, wo=True)))\n",
    "    # check if string is empty \n",
    "    if not sv_intersect_str.getvalue():\n",
    "        raise ValueError(\"No variants overlap ChIP-seq peaks\")\n",
    "    sv_intersect_df = pd.read_csv(sv_intersect_str, sep=\"\\t\", header=None).astype({3:str}).rename(columns=bed_intersect_cols)\n",
    "    if cell_type == \"all\": \n",
    "        cre_type_index = -2\n",
    "    else: \n",
    "        cre_type_index = -3\n",
    "    cre_types_list = sv_intersect_df[sv_intersect_df.columns[cre_type_index]].unique().tolist()\n",
    "    for cre_type in cre_types_list: \n",
    "        vrnt_id_intersect_cre_type = sv_intersect_df[sv_intersect_df[sv_intersect_df.columns[cre_type_index]] == cre_type].copy().vrnt_id.unique()\n",
    "        vrnt_ccre_features_df[f\"{cre_type.replace('-', '').replace(',', '')}_{cell_type}\"] = np.where(vrnt_ccre_features_df.vrnt_id.isin(vrnt_id_intersect_cre_type), 1, 0)\n",
    "    print(f\"{cell_type}: number of unique variants in dataframe: {len(vrnt_ccre_features_df.vrnt_id.unique())}\")\n",
    "# subset to CTCF features \n",
    "vrnt_ctcf_ccre_features_df = vrnt_ccre_features_df[[\"vrnt_id\", \"misexp_uniq\", \"CTCFonlyCTCFbound_all\", 'CTCFonlyCTCFbound_CD14_monocyte', 'HighCTCF_B_cell', \"HighCTCF_Neutrophil\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba7ea739-b942-4f0d-b671-571ab13519cc",
   "metadata": {},
   "source": [
    "### CpG Islands "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b00bac3e-9d00-4651-9757-3a852a599780",
   "metadata": {},
   "outputs": [],
   "source": [
    "### CpG Islands \n",
    "cpg_isl_bed = BedTool(cpg_isl_bed_path).sort()\n",
    "# identify SVs that overlap islands\n",
    "bed_intersect_cols_cpg = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\", 4:\"cpg_isl_chrom\", \n",
    "                          5:\"cpg_isl_start\", 6:\"cpg_isl_end\", 8:\"overlap\"}\n",
    "sv_intersect_cpg_isl_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(cpg_isl_bed, wo=True)))\n",
    "sv_intersect_cpg_isl_df = pd.read_csv(sv_intersect_cpg_isl_str, sep=\"\\t\", header=None).rename(columns=bed_intersect_cols_cpg).astype({\"vrnt_id\":str})\n",
    "sv_intersect_cpg_isl = sv_intersect_cpg_isl_df.vrnt_id.unique()\n",
    "# annotate variants intersecting at least one CpG island\n",
    "vrnts_cpg_island_df = all_vrnts_in_windows_df.copy() \n",
    "vrnts_cpg_island_df[\"intersect_cpg_isl\"] = np.where(vrnts_cpg_island_df.vrnt_id.isin(sv_intersect_cpg_isl), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12ed990b-a22d-4dd2-8b85-30dc3b8edbf4",
   "metadata": {},
   "source": [
    "### Chromatin States "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "94146b1d-a5d5-4990-8b4e-037024d705e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "chromhmm_15states = ['1_TssA', '2_TssAFlnk', '3_TxFlnk', '4_Tx', '5_TxWk', '6_EnhG', '7_Enh', '8_ZNF/Rpts', '9_Het', \n",
    "                     '10_TssBiv', '11_BivFlnk', '12_EnhBiv', '13_ReprPC', '14_ReprPCWk', '15_Quies']\n",
    "vrnt_bed_columns = {0:\"chrom\", 1:\"start\", 2:\"end\", 3:\"vrnt_id\"}\n",
    "bed_intersect_cols_chrom_hmm = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\", 4:\"state_chrom\", \n",
    "                                5:\"state_start\", 6:\"state_end\", 7:\"state_name\", 8:\"overlap\"}\n",
    "\n",
    "# load PBMCs chrom HMM states\n",
    "pbmcs_chromhmm_15marks_bed = BedTool(pbmcs_chromhmm_15marks_path)\n",
    "    \n",
    "# identify SVs that overlap peaks\n",
    "sv_intersect_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(pbmcs_chromhmm_15marks_bed, wo=True)))\n",
    "sv_intersect_df = pd.read_csv(sv_intersect_str, sep=\"\\t\", header=None).astype({3:str}).rename(columns=bed_intersect_cols_chrom_hmm)\n",
    "# build dataframe \n",
    "vrnts_chrom_hmm_df = all_vrnts_in_windows_df.copy()\n",
    "for state in chromhmm_15states: \n",
    "    # does state overlap variant \n",
    "    vrnt_ids_with_state = sv_intersect_df[sv_intersect_df.state_name == state].vrnt_id\n",
    "    vrnts_chrom_hmm_df[state] = np.where(vrnts_chrom_hmm_df.vrnt_id.isin(set(vrnt_ids_with_state)), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77236e1c-2094-4771-9e43-476c98ba1fd9",
   "metadata": {},
   "source": [
    "### Combine features "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "598f09c8-a6ff-44d8-a5df-fef4a84bd590",
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge different features \n",
    "dfs_to_merge = [vrnt_tad_features_df, \n",
    "                vrnt_compartment_features_df, \n",
    "                vrnt_ctcf_ccre_features_df, \n",
    "                vrnts_cpg_island_df,\n",
    "                vrnts_chrom_hmm_df]\n",
    "\n",
    "vrnt_features_merged_df = reduce(lambda  left,right: pd.merge(left,right, on=all_vrnts_in_windows_df.columns.tolist(),\n",
    "                                                              how='inner'), dfs_to_merge)\n",
    "\n",
    "# remove / in name and remove number from name for chromatin HMM columns \n",
    "rename_columns = {col:\"\".join(col.split(\"/\")).split(\"_\")[1] for col in chromhmm_15states}\n",
    "vrnt_features_merged_df = vrnt_features_merged_df.rename(columns=rename_columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d284a05f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# structural variant information \n",
    "sv_info_df =pd.read_csv(sv_info_path, sep=\"\\t\", dtype={\"plinkID\": str}).rename(columns={\"plinkID\":\"vrnt_id\"})\n",
    "\n",
    "vrnt_features_merged_info_df = pd.merge(vrnt_features_merged_df, \n",
    "                                   sv_info_df, \n",
    "                                   on=\"vrnt_id\", \n",
    "                                   how=\"left\"\n",
    "                                  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0d25c341",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n"
     ]
    }
   ],
   "source": [
    "# check for NaNs\n",
    "print(vrnt_features_merged_info_df.isnull().values.any())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "04a26b83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# results\n",
    "features_dir_path = out_dir_path.joinpath(\"features\")\n",
    "features_dir_path.mkdir(parents=True, exist_ok=True)\n",
    "# write features to file \n",
    "vrnt_features_out = features_dir_path.joinpath(\"vrnt_features_reg_annot.csv\")\n",
    "vrnt_features_merged_info_df.to_csv(vrnt_features_out, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a83fd1b3-5f29-4e5d-ba89-099184282f58",
   "metadata": {},
   "source": [
    "### Enrichment analysis \n",
    "\n",
    "* Normalise all features (including binary) and do single feature logistic regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e5b32c64-aa09-4fac-8a65-15c692bf713c",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_dict = {\n",
    "                 \"3d_genome\": ['gm12878_shared_intersect_tad_boundary', 'A_overlap', 'B_overlap', 'CTCFonlyCTCFbound_all', \n",
    "                                'CTCFonlyCTCFbound_CD14_monocyte', 'HighCTCF_B_cell', \"HighCTCF_Neutrophil\",\n",
    "                              ], \n",
    "                 \"epigenetics\":[\"intersect_cpg_isl\"] + list(rename_columns.values())\n",
    "                }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e3be963c-4927-40ee-a117-259dbf282acc",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030222\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030460\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030443\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030402\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030368\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030115\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030180\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030096\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030312\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030462\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030454\n",
      "         Iterations 107\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.029823\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030164\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030375\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030172\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030454\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030413\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030450\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030378\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030092\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030327\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030166\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030425\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049431\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047650\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048130\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049553\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048034\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048253\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047773\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046654\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.044903\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.045504\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048782\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046832\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.043372\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.045932\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.042952\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050338\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049756\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049472\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050018\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047638\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046920\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048244\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049936\n",
      "         Iterations 10\n"
     ]
    }
   ],
   "source": [
    "logistic_regr_results = {}\n",
    "logr_count = 0          \n",
    "for sv_type in [\"DEL\", \"DUP\"]:\n",
    "    input_df = vrnt_features_merged_info_df[vrnt_features_merged_info_df.SVTYPE == sv_type].copy()\n",
    "    input_df[\"svlen_norm\"] = (input_df[\"SVLEN\"] - input_df[\"SVLEN\"].mean())/input_df[\"SVLEN\"].std()\n",
    "    for key in features_dict.keys():\n",
    "        features = features_dict[key]\n",
    "        for feature in features:\n",
    "            input_df[f\"{feature}_norm\"] = (input_df[feature] - input_df[feature].mean())/input_df[feature].std()\n",
    "            y, X = dmatrices(f'misexp_uniq ~ {feature}_norm + svlen_norm', input_df, return_type = 'dataframe')\n",
    "            logit_fit = sm.Logit(endog=y, exog=X).fit(maxiter=1000)\n",
    "            log_odds, pval = logit_fit.params[1], logit_fit.pvalues[1]\n",
    "            # normal approximation confidence intervals\n",
    "            lower_conf = logit_fit.conf_int(alpha=0.05)[0][1]\n",
    "            upper_conf = logit_fit.conf_int(alpha=0.05)[1][1]\n",
    "            logistic_regr_results[logr_count] = [key, feature, sv_type, log_odds, lower_conf, upper_conf, pval]\n",
    "            logr_count += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1790e119-70d9-48f6-9eff-43e639cc3d56",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "columns_logr_enrich = [\"category\", \"feature\", \"sv_type\", \"log_odds\", \"lower\", \"upper\", \"pval\"]\n",
    "logistic_regr_results_df = pd.DataFrame.from_dict(logistic_regr_results, orient=\"index\", columns=columns_logr_enrich)\n",
    "\n",
    "# multiple testing correction\n",
    "pval_list = logistic_regr_results_df.pval.to_numpy()\n",
    "# BH FDR\n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"fdr_bh\")\n",
    "logistic_regr_results_df[\"pass_fdr_bh\"] = reject\n",
    "logistic_regr_results_df[\"pvals_corrected_fdr_bh\"] = pvals_corrected\n",
    "# Bonferroni correction \n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"bonferroni\")\n",
    "logistic_regr_results_df[\"pass_bonf\"] = reject\n",
    "logistic_regr_results_df[\"pvals_corrected_bonf\"] = pvals_corrected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a68bd17d",
   "metadata": {},
   "source": [
    "### Significant results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c793f14b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enh DUP 1.2620667347930006\n",
      "TxWk DUP 1.2483521290660107\n",
      "intersect_cpg_isl DUP 0.8790587581418782\n",
      "TssA DUP 0.8279293576899485\n",
      "TssAFlnk DUP 0.813158618928767\n",
      "ReprPC DUP 0.716611544333929\n",
      "Tx DUP 0.6961995041587603\n",
      "EnhG DUP 0.5916922200004098\n",
      "Tx DEL 0.3626878232536936\n",
      "ReprPCWk DEL 0.30791525341755954\n",
      "HighCTCF_B_cell DEL 0.2971438567503857\n",
      "TxWk DEL 0.29588000202608405\n",
      "intersect_cpg_isl DEL 0.27403255506016183\n",
      "Enh DEL 0.25040211979759047\n",
      "HighCTCF_Neutrophil DEL 0.23921334073247189\n",
      "gm12878_shared_intersect_tad_boundary DEL 0.19924272289615444\n",
      "EnhBiv DEL 0.19562396283062797\n"
     ]
    }
   ],
   "source": [
    "logistic_regr_results_pass_bonf_df = logistic_regr_results_df[logistic_regr_results_df.pass_bonf].sort_values(by=\"log_odds\", ascending=False)\n",
    "for index, row in logistic_regr_results_pass_bonf_df.iterrows(): \n",
    "    print(row[\"feature\"], row[\"sv_type\"], row[\"log_odds\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5420509",
   "metadata": {},
   "source": [
    "### Write results to file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ed966d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "category_names = {\"3d_genome\": \"3D Genome\", \"epigenetics\": \"Regulatory\"}\n",
    "\n",
    "feature_names = {'gm12878_shared_intersect_tad_boundary': \"TAD boundaries (shared)\", \n",
    "                 'CTCFonlyCTCFbound_all': \"CTCF (all)\",\n",
    "                 'CTCFonlyCTCFbound_CD14_monocyte': \"CTCF (CD14-monocytes)\", \n",
    "                 'HighCTCF_B_cell': \"CTCF (B-cells)\",\n",
    "                 'HighCTCF_Neutrophil': \"CTCF (Neutrophils)\", \n",
    "                 'intersect_cpg_isl': \"CpG islands\", \n",
    "                 'TssA': \"Active TSS\", \n",
    "                 'TssAFlnk': \"Flanking active TSS\",\n",
    "                 'TxFlnk': \"Transcr. at gene 5' and 3'\", \n",
    "                 'Tx': \"Strong transcription\", \n",
    "                 'TxWk': \"Weak transcription\", \n",
    "                 'EnhG': \"Genic enhancers\", \n",
    "                 'Enh': \"Enhancers\", \n",
    "                 'ZNFRpts': \"ZNF genes & repeats\", \n",
    "                 'Het': \"Heterochromatin\", \n",
    "                 'TssBiv': \"Bivalent/poised TSS\",\n",
    "                 'BivFlnk': \"Flanking bivalent TSS/enhancer\", \n",
    "                 'EnhBiv': \"Bivalent enhancer\", \n",
    "                 'ReprPC': \"Repressed polyComb\", \n",
    "                 'ReprPCWk': \"Weak repressed polycomb\", \n",
    "                 'Quies': \"Quiescent\", \n",
    "                 'A_overlap': \"A compartment\", \n",
    "                 'B_overlap': \"B compartment\"\n",
    "                }\n",
    "\n",
    "logistic_regr_results_df[\"feature_name\"] = logistic_regr_results_df.feature.replace(feature_names)\n",
    "logistic_regr_results_df[\"category_name\"] = logistic_regr_results_df.category.replace(category_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "7c56dda8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# all results for deletions and duplications, use Bonferroni cutoff\n",
    "logistic_regr_results_df[\"log_odds_adj\"] = np.where(logistic_regr_results_df.pass_bonf, \n",
    "                                                    logistic_regr_results_df.log_odds, \n",
    "                                                    np.nan)\n",
    "# results\n",
    "results_dir_path = out_dir_path.joinpath(\"results\")\n",
    "results_dir_path.mkdir(parents=True, exist_ok=True)\n",
    "all_results_path = results_dir_path.joinpath(\"logr_func_enrich_results_all.tsv\")\n",
    "logistic_regr_results_df.to_csv(all_results_path, sep=\"\\t\", index=False)                                      "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20e84e06",
   "metadata": {},
   "source": [
    "### Non-length adjusted enrichments "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d7e4be71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030299\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030644\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030633\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030489\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030485\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030172\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030239\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030176\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030417\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030625\n",
      "         Iterations 9\n",
      "Warning: Maximum number of iterations has been exceeded.\n",
      "         Current function value: 0.030639\n",
      "         Iterations: 1000\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.029946\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030316\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030505\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030278\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030609\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030513\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030583\n",
      "         Iterations 9\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_t/tv5/.conda/envs/tv5_base/lib/python3.7/site-packages/statsmodels/base/model.py:606: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030490\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030164\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030430\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030296\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030624\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050388\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050647\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.051551\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050535\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048839\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048971\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048330\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047158\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.045516\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046238\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.051113\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048057\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.044578\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046893\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.043429\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.052576\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050928\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050554\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.051687\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048430\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047370\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049510\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.052227\n",
      "         Iterations 10\n"
     ]
    }
   ],
   "source": [
    "logistic_regr_results_non_len_adj = {}\n",
    "logr_count = 0          \n",
    "for sv_type in [\"DEL\", \"DUP\"]:\n",
    "    input_df = vrnt_features_merged_info_df[vrnt_features_merged_info_df.SVTYPE == sv_type].copy()\n",
    "    for key in features_dict.keys():\n",
    "        features = features_dict[key]\n",
    "        for feature in features:\n",
    "            input_df[f\"{feature}_norm\"] = (input_df[feature] - input_df[feature].mean())/input_df[feature].std()\n",
    "            y, X = dmatrices(f'misexp_uniq ~ {feature}_norm', input_df, return_type = 'dataframe')\n",
    "            logit_fit = sm.Logit(endog=y, exog=X).fit(maxiter=1000)\n",
    "            log_odds, pval = logit_fit.params[1], logit_fit.pvalues[1]\n",
    "            # normal approximation confidence intervals\n",
    "            lower_conf = logit_fit.conf_int(alpha=0.05)[0][1]\n",
    "            upper_conf = logit_fit.conf_int(alpha=0.05)[1][1]\n",
    "            logistic_regr_results_non_len_adj[logr_count] = [key, feature, sv_type, log_odds, lower_conf, upper_conf, pval]\n",
    "            logr_count += 1\n",
    "            \n",
    "logistic_regr_results_non_len_adj_df = pd.DataFrame.from_dict(logistic_regr_results_non_len_adj, orient=\"index\", columns=columns_logr_enrich)\n",
    "\n",
    "# multiple testing correction\n",
    "pval_list = logistic_regr_results_non_len_adj_df.pval.to_numpy()\n",
    "# BH FDR\n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"fdr_bh\")\n",
    "logistic_regr_results_non_len_adj_df[\"pass_fdr_bh\"] = reject\n",
    "logistic_regr_results_non_len_adj_df[\"pvals_corrected_fdr_bh\"] = pvals_corrected\n",
    "# Bonferroni correction \n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"bonferroni\")\n",
    "logistic_regr_results_non_len_adj_df[\"pass_bonf\"] = reject\n",
    "logistic_regr_results_non_len_adj_df[\"pvals_corrected_bonf\"] = pvals_corrected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "58310a7d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enh DUP 1.3483934099463208\n",
      "TxWk DUP 1.3371875314903858\n",
      "intersect_cpg_isl DUP 0.9900244119091234\n",
      "TssA DUP 0.9199302894668124\n",
      "TssAFlnk DUP 0.914645477756096\n",
      "HighCTCF_B_cell DUP 0.9067217923366574\n",
      "ReprPC DUP 0.818066833360135\n",
      "Tx DUP 0.8038045067062993\n",
      "HighCTCF_Neutrophil DUP 0.7961936907052427\n",
      "CTCFonlyCTCFbound_CD14_monocyte DUP 0.719790682335864\n",
      "EnhG DUP 0.6773341661297582\n",
      "EnhBiv DUP 0.6090085756204552\n",
      "gm12878_shared_intersect_tad_boundary DUP 0.521749321061759\n",
      "TssBiv DUP 0.470566765240087\n",
      "Tx DEL 0.3771288931606897\n",
      "HighCTCF_B_cell DEL 0.329535373226274\n",
      "ReprPCWk DEL 0.32921977095901356\n",
      "TxWk DEL 0.30888578279639145\n",
      "intersect_cpg_isl DEL 0.30070858471699463\n",
      "Enh DEL 0.27679842065477195\n",
      "TxFlnk DUP 0.2751730516402117\n",
      "HighCTCF_Neutrophil DEL 0.27019044147293353\n",
      "gm12878_shared_intersect_tad_boundary DEL 0.2269056914946834\n",
      "EnhBiv DEL 0.21838755116601496\n",
      "TssA DEL 0.20067859853742478\n"
     ]
    }
   ],
   "source": [
    "logistic_regr_results_non_len_adj_pass_bonf_df = logistic_regr_results_non_len_adj_df[logistic_regr_results_non_len_adj_df.pass_bonf].sort_values(by=\"log_odds\", ascending=False)\n",
    "for index, row in logistic_regr_results_non_len_adj_pass_bonf_df.iterrows(): \n",
    "    print(row[\"feature\"], row[\"sv_type\"], row[\"log_odds\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a12a8e82",
   "metadata": {},
   "outputs": [],
   "source": [
    "logistic_regr_results_non_len_adj_df[\"feature_name\"] = logistic_regr_results_non_len_adj_df.feature.replace(feature_names)\n",
    "logistic_regr_results_non_len_adj_df[\"category_name\"] = logistic_regr_results_non_len_adj_df.category.replace(category_names)\n",
    "\n",
    "# all results for deletions and duplications, use Bonferroni cutoff\n",
    "logistic_regr_results_non_len_adj_df[\"log_odds_adj\"] = np.where(logistic_regr_results_non_len_adj_df.pass_bonf, \n",
    "                                                                logistic_regr_results_non_len_adj_df.log_odds, \n",
    "                                                                np.nan)\n",
    "all_results_no_len_adj_path = results_dir_path.joinpath(\"logr_func_enrich_results_all_no_len_adj.tsv\")\n",
    "logistic_regr_results_non_len_adj_df.to_csv(all_results_no_len_adj_path, sep=\"\\t\", index=False)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f4872a3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
